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Abstract. We report a new theory of dissipative forces acting between colliding viscoelastic bodies. The 
impact velocity is assumed not to be large, to avoid plastic deformations and fragmentation at the impact. 

The bodies may be of an arbitrary convex shape and of different materials. We develop a mathematically 
rigorous perturbation scheme to solve the continuum mechanics equation that deals with both displacement 
and displacement rate fields and accounts for the dissipation in the bulk of the material. The perturbative 
solution of this equation allows to go beyond the previously used quasi-static approximation and obtain the 
dissipative force. This force does not suffer from the physical inconsistencies of the latter approximation 
and depends on particle deformation and deformation rate. 

PACS. 45.50.Tn Collisions - 45.70.-n Granular systems - 46.35.-l-z Viscoelasticity, plasticity, viscoplas¬ 
ticity 


1 Introduction 

Granular materials are abundant in nature; they range 
from sand and powders on Earth to planetary rings and 
dust clouds in outer space jl|2|3|4 |5|. These material ex¬ 
hibit very unusual properties, demonstrating solid-like, 
liquid-like or gas-like mm behavior depending on the 
external load or magnitude of agitation [10111112] , The 
physical reason for many unusual phenomena in granular 
media is the nature of inter-particles interactions in these 
systems. Contrary to molecular or atomic systems, where 
particles interact only trough conservative, elastic forces, 
the interaction between granular particles include dissi¬ 
pative forces. This happens because the grains are them¬ 
selves macroscopic bodies, which contain macroscopically 
large number of microscopic degrees of freedom. Hence, 
during an impact of such bodies their mechanical energy, 
associated with the translational or rotational motion, or 
with the elastic deformation of the particles, is partly 
transformed into the internal degrees of freedom, that is, 
into heat. In many applications however, the temperature 
increase of the grains is insignificant and may be neglected 
[6]. Obviously, for an adequate description of granular me¬ 
dia one needs a quantitative model of inter-particles forces, 
which includes both elastic and dissipative interactions. 

The elastic part of the inter-particle force is known for 
more than a century from the famous work of Hetrz [13] . 
Hertz obtained a mathematically rigorous result for the 
force acting between elastic bodies at a contact, provided 
the deformation of the bodies is small as compared to 


their size; the theory has been developed for particles 
of an arbitrary convex shape. In spite of a large impor¬ 
tance for applications, the rigorous derivation of the dis¬ 
sipative force is still lacking. Presently there exist phe¬ 
nomenological expressions for the dissipative force, which 
exploit either linear, e.g. [14115] or quadratic [TS] depen¬ 
dence of the force on the deformation rate. Neither lin¬ 
ear, nor quadratic dependence, however, is consistent with 
the experimental data, e.g. mu- A derivation of the 
dissipative force from the first-principles has been under¬ 
taken in Ref. [181 . A very restrictive approximation used 
in this work - the assumption that only the shear defor¬ 
mation is important, substantially undermines its appli¬ 
cation. A complete derivation of the dissipative force be¬ 
tween viscoelastic bodies from the continuum mechanics 
equations has been done only recently [19] within a quasi¬ 
static approximation. In this approximation it is assumed 
that the displacement field in the bulk of colliding bod¬ 
ies completely coincides with that for a static contact [T9i • 
The correct functional dependence of the dissipative force, 
derived in Ref. [T9] has been already suggested (without 
any rigorous mathematic analysis) in the earlier work of 
Kuwabara and Kono E03- In the later studies HTTHl a 
flaw in the derivation of the dissipative force of Ref. [T9] 
has been corrected. Still the restrictive assumption of the 
quasi-static approximation has been exploited mrm . 

Physically, the quasi-static approximation assumes the 
immediate response of the particle’ material to the exter¬ 
nal load. Two conditions are to be fulfilled in order to 
make this approximation valid: (i) the characteristic de- 
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formation rate should be much smaller than the speed of 
sound in the system and (ii) microscopic relaxation time 
of the particle’s material should be much shorter than the 
duration of the impact. The microscopic relaxation time 
quantifies the time needed for the material of a deformed 
body to respond to the applied load; in what follows we 
will give the detailed definition of this quantity. To go be¬ 
yond the quasi-static approximation, that is, to take into 
account the deviation of the displacement field in the bulk 
of a deformed body from the static displacement field, we 
develop a perturbation approach based on small parame¬ 
ter — the ratio of microscopic relaxation time and collision 
duration. In the most of important applications this ra¬ 
tio is indeed small, which implies that for the first time 
we rigorously derive a dissipative force acting between 
viscoelastic particles. Although the quasi-static approxi¬ 
mation is based on the physically plausible approach, it 
possesses some inconsistency. This inconsistency is not so 
visible for a collision of particles of the same material. At 
the same time when particles of different materials suffer 
an impact, the quasi-static approximation predicts non¬ 
equal dissipative forces acting between the bodies, which 
definitely violates the third Newton’s law. The other in¬ 
consistency is related to the dependence of the dissipative 
force on the Poisson ratio - within the quasi-static ap¬ 
proximation one obtains zero dissipative force for the case 
of vanishingly small elastic shear module; this is definitely 
not physical. These difficulties of the quasi-static approx¬ 
imation are discussed in detail below. 

Our new theory, based on the perturbation scheme, 
is mathematically rigorous and the obtained dissipative 
force is free from the above inconsistencies. While in the 
present work we analyze a general case of an impact of 
viscoelastic bodies of an arbitrary shape and of different 
materials, the results for a more simple case of a collision 
of a sphere with un-deformable plane has been reported 
earlier [23] . 

The rest of the paper is organized as follows. In the 
next Sec. II we introduce the equation of motion of vis¬ 
coelastic medium which we solve for the case of interest 
in the next sections. In Sec. Ill the solution for the static 
contact is considered; here we illustrate the general ap¬ 
proach and derive the classical Hertz law. In Sec. IV the 
dynamic contact is addressed. We elaborate the perturba¬ 
tion scheme and using this scheme derive in Sec. V the 
next-order solution. In Sec. VI we present our new the¬ 
ory for the dissipative force between colliding viscoelastic 
bodies and finally in Sec. VII we summarize our findings. 


2 Equation of motion for viscoelastic medium 


governed by the equation for a continuum medium which 
reads, e.g. [24;, 

pu = V ■ a = V • ( a el + a v ) . (1) 

Here p is the material density, u = u(r) is the displace¬ 
ment field in a point r and a is the stress tensor, comprised 
of the elastic a el and viscous part a v . The elastic stress 
linearly depends on the strain tensor, 

Uij = i (Vj uj + V jut ), 

built on the displacement field [24]: 


<J(u) = 2 E 1 




-j- E^SijUll . 


( 2 ) 


Similarly, the viscous stress linearly depends on the strain 
rate tensor [21] : 


o«(u) = 2r?i 




+ V2^ijUll ■ 


(3) 


Here E\ = 2 (i+v) an< ^ ^ 2 = 3 (i- 2 „) i with Y and v being 
respectively the Young modulus and Poisson ratio of the 
body material. ?yi and rj 2 are the viscosity coefficients for 
the shear and bulk viscosity and i,j,l denote Cartesian 
coordinates; the Einstein’s summation rule is applied. 

The elastic deformation implies that, after separation 
of the contacting particles, they completely recover their 
initial shape so that no plastic deformation remains. Only 
such deformations will be addressed below. 


3 Static contact. Hertz theory. 

To introduce the notations and illustrate the main tech¬ 
nical ideas, we start from the simplest case of a static 
contact, that is, we consider a time-independent contact 
of two convex bodies. We assume that only normal forces, 
with respect to the contact area, act between the parti¬ 
cles. We place the coordinate system in the center of the 
contact region, where x = y = z = 0 (Fig. [T]). Let the 
displacement field in the upper body, located at z > 0, 
be u(r), while in the lower body, located at z < 0 be 
w(r). Then the deformation £ which is equal to the sum 
of the compressions of the both bodies in the center of 
the contact zone is related to the z-components of the 
displacements of the upper and lower bodies’ surfaces at 
the contact plane u z (x, y, 0) and w z (x, y , 0), Fig. [1] It may 
be shown m that for the bodies of arbitrary shape the 
following relation holds true: 


When two viscoelastic bodies are brought in a contact, 
so that they are deformed, an interaction force between 
the bodies arise. Generally, it contains elastic and viscous 
parts; for a static contact however, only the elastic force 
appears. To compute the forces, one needs to find a stress 
that emerges in the bodies and integrate the stress over the 
contact zone. The distribution of stress in the material is 


Bix 2 + B 2 y 2 + u z {x , y, 0) + w z (x , y, 0) = £, (4) 


where the constants B\ and B 2 are related to the radii of 
curvature of the bodies’ surfaces near the contact [24] . 


2 ( B i + B 2 ) - — + — + — + — , 


Rl I?2 


(5) 
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Fig. 1. Illustrates a simple case of a collision of two visco¬ 
elastic spheres (the dashed profiles show undeformed bodies) 
in the according coordinate frame. Note however, that in the 
text a general case of arbitrary convex bodies is addressed. 



Here i?i, R 2 and R .\, R' 2 are respectively the principal 
radii of curvature of the first and the second body at the 
point of contact and ip is the angle between the planes 
corresponding to the curvature radii R\ and R[. Equa¬ 
tions and © describe the general case of the con¬ 
tact between two smooth bodies (see [24] for details). The 
physical meaning of ® is easy to see for the case of a 
contact of a soft sphere of a radius R (R± = R 2 = R) with 
a hard, undeformed plane (R.\ = R' 2 = 00 ). In this case 
B\ = B 2 = 1/21?, the compressions of the sphere and of 
the plane are respectively u 2 (0,0, 0) = £ and w z = 0, and 
the surface of the sphere before the deformation is given 
by z(x,y) = ( x 2 + y 2 )/2R for small 2 . Then ([3]) reads in 
the flattened area u z (x,y) = £ — z(x,y ), that is, it gives 
the condition for a point z(x, y) on the body’s surface to 
touch the plane z = 0. While Eq. a defines the displace¬ 
ment on the contact surface, the displacement fields in 
the bulk of the first (upper) and second (lower) bodies are 
determined by the following equations. 

V ■ a el (u) = 0, V ■ a el (w) = 0. (7) 

Both equations may be solved by the same approach, 
therefore in what follows we consider the solution for the 
upper body with z > 0. Using Eq. © which relates the 


stress and strain tensors, we write: 


Vj-ofj — E\Aui + ( E 2 + 77 E\ ] V -i V/ Uj — 0 


( 8 ) 


where the elastic constants refer to the upper body (for 
the notation simplicity we do not add now the additional 
index specifying the body - it will be done later). 

To solve the above equation we use the approach of [33] 
and write the solution as 

U (0) = / (0) e 2 + , (9) 

where = K^zf^ +ip(°\ K is some constant to be 
found and and ip^ are unknown harmonic functions. 
We assume the lack of tangential stress at the interface, 
which is e.g. fulfilled when the bodies at a contact are of 
the same material. Taking into account that 

Au = AV<pW = 2K^\7^E- (10) 

dz 

and 

Vu = (1 + 2A'(°))^^, (11) 

OZ 

we recast Eq. ([5]) into the following form: 




2E 1 K (0) + 

(1 + 2/i (0) ) ( 


( 12 ) 


E 2 + 


Ei 


0 /(°) 

1 dz 


= 0, 


which implies that 

K (o) = 1 3E 2 + E! ' 

2 3 E 2 + 4 E\ 


(13) 


Consider now the boundary condition for the stress tensor. 
Obviously, on the free boundary all components of the 
stress vanish. In the contact region, located at the surface, 
z = 0, the tangential components of the stress tensor af x 
and <rf v vanish as well, while the normal component of the 
stress tensor reads, 

n-a el = -af z = P z (14) 


where n = (0, 0, —1) is the external normal to the upper 
body on the contact plane and P z is the normal pressure 
acting on the contact surface. Therefore the boundary con¬ 
ditions have the following form: 


L=o 


el I 

zy lz=o 


= 0 ; 


12=0 


= -P,. 


(15) 


Using the expression © for the elastic part of the stress 
tensor, together with the displacement vector © we recast 
the boundary conditions m into the form: 


d_ 

dx 

d_ 

dy 

d_ 

dz 


3 ^i f (o) . 

4£i + 3 E 2 dz 

3 ^i f (o) | ofrfr 

AEi + 3E 2 dz 


3Ei 


AE\ + 3E 2 


/( 0 , + 2 |) 


= 0 

2=0 

= 0 

2=0 

Pz 


(16) 

(17) 

(18) 
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From equations (USD and m follows the relation between 
/(°) and §§ at 0 = 0: 


d4> , 3 Ei 

~dz + 2 4Ei + 3E 2 


/ (0) 


= const = 0 . 


(19) 


z =0 


dijj 

dz 


E, 


-/ (0) ■ 


2 -b 3E 2 
Substituting the last relation into (flSll yields 


df(o) 


dz 


z=0 


4 1'J'i + 3i?2 

ifi (Ei 4 - 3E2) 


P 2 . 


( 21 ) 


/ (0) (r) = 


4_£7i -b 3E 2 
‘ZttE\(Ei -b 3E 2 } 


P z {x',y') dx'dy' 


r - r' 


( 22 ) 


■°zlz=0 


=0 = (1 + ^ (0) ) / 


( 0 ) 


dip 

+ 

z =0 02: 


z=0 


which together with (BUI) and definition of K^°\ Eq. (1131) 
yields, 

1 


= - f (0) 
\z=0 o J 


z =0 


(23) 


u z {x,y,z = 0) = -- 


(l-^i 2 ) ff af z {x',y',z = d) dx'dy’ 


nYi 


w z (x, y,z = 0) = - 


(! ~ v l 

nY 2 


Hence, with Eq. (fTUl) the relation (HD takes the form: 

1 fl-v 2 , , ff P z (x',y') 

Y 2 


Vi 


-dx'dy’ = 

|i- - 1 - | 

= £ ~ B lX 2 - B 2 y 2 , (26) 


The constant in the above relation equals to zero, since it 
holds true independently on the coordinate that is, also at 
the infinity; at the infinity, however, the deformation and 
thus the above functions vanish. Since f(°\ ip as well as 
dip/dz are the harmonic functions, the condition that their 
linear combination vanishes on the boundary, Eq. (USD , 
implies that it is zero in the total domain, that is, 


The equation (1261) is an integral equation for the unknown 
function P z (x,y). We compare this equation with the math¬ 
ematical identity m 


dx'dy' 
S Ir-r'l 


( 20 ) = ^ 

2 


x' 2 y' 2 

Jl 


(27) 


0 


1 -- 


X 


dt 


a 2 +t b 2 +t J ^/(a 2 + t)(b 2 + t)t ’ 


Since /is a harmonic function, one can use the relation 
between the normal derivative of a harmonic function on 
the surface and its value in the bulk, as it follows from the 
theory of harmonic functions (see e.g. mm, hence we 
find: 


where integration is performed over the elliptical area 
x' 2 /a 2 + y' 2 /b 2 < 1. The left-hand sides of the both equa¬ 
tions, m and ([SHD , contain integrals of the same type, 
while the right-hand sides contain quadratic forms of the 
same type. Therefore, the contact area is an ellipse with 
the semi-axes a and b and the pressure is of the form 


where S is the contact area. Using Eq. © we can write 
^-component of the zero-order displacement at 0 = 0 as 


P z (x,y) = const \J\ — x 2 /a 2 — y 2 /b 2 . 

The constant here may be found from the total elastic 
force E e i acting between the bodies. Integrating P z (x,y) 
over the contact area we get E e i, which then yields the 
constant. Hence we obtain 


Pz{x,y) = 


3E et 
27 Tab 


x 2 y 2 


b 2 ‘ 


(28) 


We substitute (B51) into (BUI) and replace the double in¬ 
tegration over the contact area by integration over the 
variable t, according to the above identity. Thus, we ob¬ 
tain an equation containing terms proportional to x , y 2 
and a constant. Equating the corresponding coefficients 
we obtain 


If we now express E\ and E 2 in terms of v\ and Yi, where 
and Yi are the according constants for the upper body 
(recall that we consider the upper contacting body) one 
obtains from Eqs. B3l) . B2l) and m-- 


e= 


El = 


E eI E r°° dt 

n Jo \J (a 2 + t)(b 2 + t)t 

F e \D f°° dt 


FeX N( C) 


(29) 


7T Jo (a 2 + t)yj (a 2 + t)(b 2 + t)t 
F e \D M{ C) 


( 24 ) 

The same considerations may be performed for the lower 
body. Taking into account that the external normals for 
the upper and lower bodies as well as the exerted pressures 
are equal up to a minus sign (n up = — nj, P z = P z<up = 
—P z ,\), we obtain, 


E 2 = 


7r a 2 b 
F e \D f°° 


dt 


7T Jo {b 2 + t)\J (a 2 + t){b 2 + t)t 
= FeiD M{ I/O 
7 r ab 2 ’ 


(30) 


(31) 


af z (x\ y', 0 = 0) dx'dy' where 


3 / 1 -^ 

4 l Y 1 


1-^2 

y 2 


(25) 


(32) 
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and £ = a 2 /b 2 is the ratio of the contact ellipse semi-axes. In this case Co = 1, N( 1) = 7r, and M( 1) = n/2, leading 
In (151) (15T1) we introduce the short-hand notation^] to the solution of (1571) . (1581) : 


N(0 = 


M( C) = 



d t 

\/(i + (t)(i+t)t 
d t 

(1 + t)\/(l + t)( 1 + £ t)t 


(33) 

(34) 


From the above relations follow the size of the contact 
area, a, b and the deformation £ as functions of the elastic 
force F e \ and (known) geometric coefficients B\ and S 2 . 

The dependence of the force F e i on the deformation 
£ may be obtained from scaling arguments. If we rescale 
a 2 —y era 2 , b 2 ab 2 , £ —> a£ and F e i —> a 3 ^ 2 F e \, with a 
constant, Eqs. (HE) (ED) remain unchanged. That is, when 
£ changes by the factor a, the semi-axis a and b change 
by the factor a 1 / 2 and the force by the factor a 3 / 2 , i.e., 
a ~ £ 3 / 2 , b ~ £ 3 / 2 and 

F e i = const £ 3 ^ 2 . (35) 


The dependence (1551) holds true for all smooth convex bod¬ 
ies in contact. To find the constant in (l35l) we divide (1511) 
by (1501) and obtain the transcendental equation 


b 2 vW 1/0 

5T = ~m~ (36) 

for the ratio of semi-axes £. Let £o be the root of Eq. ED, 
then a 2 = £o b 2 and we obtain from Eqs. ED, ED: 


FelDNj Co) 

7r b 


Bx 


F e \D M(Co) 

7T £ 0 6 3 


(37) 

(38) 


where iV(£o) and M(£o) are pure numbers. Equations (1571) . 
(1551) allow us to find the semi-axes b and the elastic force 
F e \ as functions of the compression £. Hence we obtain the 
force, that is, we get the according constant in Eq. ED 

m-- 






(39) 


Similarly we can relate the deformation £ and the semi¬ 
axes a of the contact ellipse [55] : 


a = 


( M(Co) V /2 1/2 

Kn^bJ ? 


(40) 


Note that £o is a constant determined by the collision ge¬ 
ometry. 

For the special case of contacting spheres of the same 
material (a = b), the constants B\ and B 2 read 


< 41 » 

1 The function N(Q and M(£) may be expressed as a com¬ 
bination of the Jacobian elliptic functions _E(£) and R'(C) [28 . 


a 2 = R eS £ (42) 

F el = P e /2 ; p = 3 (1 2 ^ 2) ^ ; (43) 

where we use the definition (1551) of the constant D. This 
contact problem was solved by Heinrich Hertz in 1882 [ID- 
it describes the force between elastic particles. For inclas- 
tically deforming particles it describes the repulsive force 
in the static case. 


4 Dynamical contact. Perturbation scheme 

For the most important applications the viscous forces are 
significantly smaller than the elastic forces and the bod¬ 
ies material is rigid enough to neglect inertial effects for 
collisions with not very large velocities. Let us estimate 
the magnitude of different terms in Eq. ©■ This may be 
easily done using the dimensionless units. For the length 
scale we take R, which corresponds to the characteristic 
size of colliding bodies, while for the time scale we use r c 
- the collision duration. Then vo = R/t c is the charac¬ 
teristic velocity at the impact. Taking into account that 
differentiation with respect to a coordinate yields for di¬ 
mensionless quantities the factor 1/i?, and with respect to 
time - 1/t c , we obtain 

- A! X7a el Ai = T rel /r c , (44) 

pu ~ pw ~ A 2 Va el A 2 =Vq/c 2 . (45) 

Here c 2 = Y/p and T re i = rj/Y characterize respectively 
the speed of sound and the microscopic relaxation time in 
the material and p ~ 771 ~ 7? 2 [T9]. 

Neglecting terms, of the order of Ai and A 2 we get 

V • <r e '(u) = 0, V • cr el (w) = 0, (46) 

which yields the static displacement fields u = u(r) and 
w = w(r). This approximation corresponds to the quasi¬ 
static approximation, used in the literature I19I21I22I25I26I . 
Neglecting terms of the order A 2 but keeping these of the 
order of Ai, leads to the following equation 

V ■ (a el {u) + d B (u)) =0, V • (a el { w) + ^(w)) = 0. 

(47) 

That is, to go beyond the quasi-static approximation one 
needs to find the solution of Eq. ED which contains both 
the displacement fields u, w, as well as its time derivatives, 
u, w. Eq. (1571) needs to be supplemented by the boundary 
conditions. These correspond to vanishing stress on the 
free surface and given displacement in the contact area. 

In a vast majority of applications Ai = T re i/T c <C 1, 
which implies that the viscous stress is small as compared 
to the elastic stress. This allows to solve Eq. (1571) pertur- 
batively, as a series in a small parameter. Here we follow 
the standard perturbation scheme, e.g. [5j: To notify the 
order of different terms we introduce a ” technical” small 
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parameter A, which at the end of computations is to be 
taken as one. Hence one can write, 


With this notations and using Eqs. 0, 0 and m 
write: 


we 


a = d (0 ) + A<x < ' 1) + A 2 a < ' 2 ' 1 + ... 
and respectively, 

u(r) = u (0) (r) + Au« (r) + A 2 u (2) (r) + ..., 
w(r) = (r) + Aw^ (r) + A 2 w^ 2) (r) + ... 


(48) 

(49) 

(50) 


(“ <0) ) 


= a 


'(i) 


Vl -el( 0) 
Ei 


(54) 


+ v 2-m 


Ei 

El 


(1 + 2 K^)^-Si. 


and accordingly the divergence of this tensor: 


Substituting the expansions BSD and BHD, BHD into 
Eq. BID yields a set of equations for different order in A. 
Zero-order equations with the according boundary condi¬ 
tions read, 




2r ?1 ^°) + (l + 2^ 0 ))(r ?2 + |-) 


V, 


= 0 . 


B lX 2 


( U (0) ) 

B 2 y 2 + u^\x, y,0) 


(w (0) )=0, (51) 


3(E ir i2-E 2m ) dfW 

" V i - 


(4Ei + 3 E 2 ) 


dz 


df(o) 

dz 

(55) 


«4 0) (+?/, 0) = £, 


while the first-order equations with the boundary condi¬ 
tions have the form 


( 


•el ( 


cr“"(u^ 
a e ‘(w (1) 


= 0. 


4 1} (+ 


y, o) 


)++(u(°) 

) + ff’(w (0 > 

- fflW (a;, y, 0) = 0 


where Eqs. BSD, m and Eq. USD for have been 

used. If we now apply Eq. (1211) for df^ /dz and again 
Eq. (fT51) for the constant K^ 0 \ we find the 22 -component 
of the first-order dissipative tensor on the contact plane, 
2 = 0: 


= 0, 


(52) 


a v z ^\x,y,0) = aa e z l z io) (x,y,0) 

3rj 2 + rji 

a =-. 

Ei + 3E 2 


(56) 

(57) 


and so on. Note that the zero-order equation m cor 
responds to case of a static contact which has been con¬ 
sidered in detail above. This also corresponds to the quasi¬ 
static approximation widely used in the literature, e.g. [191211221251261 . 
Also note that in the proposed perturbation scheme, only <+ = 

zero-order problem (1511) has non-zero boundary conditions, 
corresponding to the boundary conditions 0 of the ini¬ 
tial problem; all other, high-order perturbation equations, 
have homogeneous boundary conditions. Such partition of 
the boundary conditions is justified due to the linearity of 
the problem. 

Note that for the zero-order solution the condition 
erf(,(u (C 9) = er e, (w(°)) is fulfilled at the contact plane 2 = 

0, as it directly follows from the construction of the so¬ 
lution. For the first-order solution, however, we need to 
additionally request the condition for the first-order stress 
tensor: 


Similar relation may be obtained for the lower body. Using 
the definitions of E\ and E 2 the coefficient a reads for each 
of the bodies, 


(l + t/0(l-2i/Q 
Yi 


2% (i) + 2 ? 7l(i) 


(58) 


< z (u ,0) )+4( u 


(i) 


2—0 


(<Jw (0) )+4(w (1} )) z=o , (53) 


which implies the equivalence of the first-order stress ten¬ 
sor, expressed in terms of deformation and deformation 
rate of the upper body and of the lower one. 


where the subscript i = 1,2 specifies the body - i = \ for 
the upper body and i = 2 for the lower one. The above 
relation corresponds to the according approximation of 
Ref. [211122] and coincides with the result of [21122] , where 
the necessary corrections have been introduced. Note, how¬ 
ever, that quasi-static approximation occurs to be incon¬ 
sistent for the case of contact of particles of different ma¬ 
terial: Indeed, the condition BHD is possible only if the 
first-order elastic terms are taken into account. Obviously, 
this may not be achieved within the quasi-static approx¬ 
imation, which uses only the first-order dissipative stress 
ali 1 ' 1 . The values of on the contact plane are different 
for the upper and lower body for different materials, since 
op ^ a 2 [see Eqs. BHD-BHD]: that is, the third Newton’s 
law for this case is violated. 

Consider now the first-order equation B2D: 


VjO 


el( 1 ) 


+ < 4 f) = o. 


(59) 


5 First-order solution. Beyond quasi-static 
approximation. 

Again we will consider the upper body with 2 > 0 and 
introduce, for convenience, the following notations: 



Due to the linearity of the problem, one can represent 
the first-order displacement field as a sum of two parts, 
uW = u (1) + uW, which correspond to the two parts of 

the elastic tensor, ofE = ofE (u^) + dfE (u^). Here 

the first part of is the solution of the inhomogeneous 
equation with homogeneous boundary conditions: 


v,^ (1) = 


- T el{ 1) 


z=0 


= a el ^ 
yz 


2=0 


= a eZ(1) 

^ zz 


= 0. 


2=o 


(60) 

(61) 
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while the second part is the solution of the homogeneous 
equation, 

(62) 


V,^ (1) =0, 


and the third boundary condition, dfS 1 '* = 0 at z = 0 is 
automatically fulfilled. Taking into account that function 
/(D 

is harmonic, we obtain, 


with a given first-order displacement u z at the contact 
plane; this is to be obtained from the boundary condition 
(1551) and consistency condition (1551) . The boundary prob¬ 
lem (1621) is exactly the same as the above problem m 
for the zero-order functions. Hence the same relation 
holds true for the first-order functions, that is, 


V,-c£? (1) = -£fiV,-^ (1) 


' 3 u ij 


dz 


vP 


(1 ~ v i) 
z=0 7 tYi 

dfP(u^(x',y',z = 0))dx l dy 1 


(63) 


Using the above equation together with Eq. dSSD we recast 
Eq. (15U1) into the form, 

df^__ 3(^277! -E, rj 2 ) ^ a/<°) 

1 1 dz (4Ei+3E 2 ) 1 dz 

which implies the relation between functions f 1 ' 1 ' 1 and 


To solve Eq. (1501) we write the displacement field u^ 1 ) 
in a form, similar to this of the zero-order solution ([5jl : 


u (1) = f {1) e z + V(/3 (1) 


(64) 


where <p W = K^zf^ +ip( 1 \ is some constant and 
/W and are harmonic functions. Then we can write 


the stress tensor a. 


■ 2 ( 1 ) 


as 


/(i) = -/3/( o) 

_ 3 (E 2 m - Ei 772 ) 
E 1 (3E 2 +4E 1 ) ' 

Using Eq. (1551) with =-|we write for 
z 2 J 2 dz ’ 

substituting there from Eq. (1711) we arrive at 


dfl 1] = (1 + 2A«) +J«Vj/W) + 


4 1} = -r/ 9 f/ (0) - 


a/W' 

92 : 


(71) 

(72) 


(73) 


(74) 


_£/2 — —.£/i 

3 


9/ (1) 

a* ^ 


2 EiV,V^ (1) . 


+ 2EiA' (1) zV,V J -/ (1) 
(65) 


If we choose K'E) = — 1 the above stress tensor takes the 
form 


where /(°) is given by Eq. (1551) . Thus, the above relation 
presents the solution for the displacement ui Taking 
now into account the relation (1551) between and 'u| 0) 
at the contact plane, as well as the expression (1241) for u 
there, we find for at z = 0: 


(0) 


o-f! (1) = + 2S 1 V i V j V’ (1) (66) -( 1 ) ( X ~ ^ 1 ) ff ffi d e z l z {0) (x f , y’, z = 0) dx'dy' 


7rFl 


r - r' 


(75) 


and the boundary conditions ED read: 


ei(l) 




_ d_ 

:=o <9x \ dz y 


= 0 , 


z=0 


2=0 


d / ch/^ 1 ) \ 
dy \ dz ) 


= 0 . 


z=0 


Therefore we conclude, 

<9'0W 


92: 


= const = 0 , 


(67) 

( 68 ) 

(69) 


where the subscript ”1” indicates that the constant /3\ 
refers to the upper body. Similar considerations may be 
done for the lower body, z < 0, yielding: 




z=0 


nY 2 

1} (w^ 1 ) (x',y',z = 0)) dx'dy' 


(76) 


and 


z=0 


where the last equation follows from the condition that 
ip^ vanishes at the infinity, x, y —> 00 , where the defor¬ 
mation is zero. Since is a harmonic function, we con¬ 
clude that the vanishing normal derivative on a boundary, 
Eq. (1691) . implies that the function vanishes everywhere, 
that is, ip'- 1 \x,y, z) = 0 (see e.g. ED- Hence 


^(1) = (1 ~ UP ['f' P2crt l z 0 ) (x',y',z = 0 ) dx'dy' ^ 


7 tY 2 


r — r' 


Now we apply the consistency condition (1551) , using Eq. (1551) 
for the both bodies, 


2=0 


a£ (1) = -EizViVjfM 


(70) 


(aiif+af ) (u< 1 >)) 

= {a 2 a e zz ^ + a-fl (1) (w (1) )j 


z=0 


(78) 
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where we also take into account that the following parts 
of the stress tensor vanish on the contact plane: 


Eq. (1751) then yields, 


2=0 


= ^ e i ( 1 ) (w (1) ) 


= o. 


2=0 


^ ( 1 ) (wW) 


z=0 


— (oq — 0:2) cr 


3/(0) 


2=0 


+ ^ ( 1 ) (Q (1) ) 


2 = 0 


(79) 


Now we use the boundary condition ED, 

+ W™ = «(!) + fiW + ^(1) + tfU) = 0 , 

and applying Eqs. (1M1) . (1751) . (1571) and (1771) for uP, zl£ 1 ' > , 
to* 1 ) and wP we obtain, 


(JhDi + fcD 2 )afW _ Dl afP( u«) 


-D 2 a e z f\ w (1) )) 


dx'dy' 


* = °> 


z=o |r — r , 

where we introduce the short-hand notations, 

A = (1 - vf)/Yi, * = 1,2. 

From the above equation, together with Eq. ed follows 
the relation for the first-order elastic tensor: 


afy\ u«) = ( ^i g i +^2 

22 V 2=0 V A + D 2 
D 2 {u\ — 0 : 2 ) 
A + -A 


(80) 


e2(0) 


2=0 


Finally we obtain, taking into account that the total first- 
order stress on the contact plane is a sum of two parts - 
the elastic one, given by Eq. ED- and the dissipative part 
from Eq. (1551) . 


T (i) 


where 


2=0 


= KP 


T ei(l) 


2 = 0 


= A & el j 0) 


2=0 


A = 


(ai + /3i)A + (a 2 + [3 2 )D 2 


(81) 


(82) 


A + D 2 

Again we take into account that the component &f!P (u fl )) 
of the stress tensor vanishes on the contact plane. The con¬ 
stant A may be written, using Eq. (1551) and m for ao(i) 
and ai(j) as 

7i A + 72 D 2 


A = 


7 i = 


A + D 2 


1 + 


(83) 


Yi 


3 » 7 i(i)(l A v?) + %(<)(! - 2 Vif 


Using the above Eqs. (1751) . (1571) and ED we can write 
the explicit expression for the viscous pressure Pp 1 ' 1 = 
acting between the colliding bodies: 


^ o P 


2=0 


P^\x,y) = - 


3A 


*DN(( 0 ) ^a2-(x 2 +y 2 Co) ’ 


(84) 


where a depends on £ according to Eq. (SOD and all other 
notations have been introduced in the previous section. 


6 Dissipative Force 

Now we can write the dissipative force acting between 
particles. It corresponds to the force associated with the 
viscous constants, that is, with the first-order stress ten¬ 
sor aP . Integrating this stress over the contact area, we 
obtain, 

pv(i) = e-CD (a;, y)\ z=0 dxdy , 

so that Eq. (1511) yields, 

A (1) = -AF 2 ei(0) , (85) 

where is the normal force corresponding to the elas¬ 

tic reaction of the medium. It is equal to the Hertzian 
force, Eq. ED; taking the time derivative of this force we 
finally obtain: 


FZ<»=-lACoVs£. ( 86 ) 

Here the constant Co, defined by Eq. ED- is determined 
by the geometry of the colliding bodies and their material 
properties (see the discussion after Eq. ED)- 

Hence the total force acting between two viscoelastic 
bodies reads in the linear approximation with respect to 
the dissipative constants: 

Ftot = Co? /2 - \ACoVU , (87) 

where the relation between the deformation £ and the axis 
a of the contact ellipse is given by Eq. ED as in the static 
Hertz theory. Note however, that contrary to the Hertz 
theory the size of the contact ellipse is determined now not 
by the total force acting between the bodies, but by the 
elastic part of the total force, Aot + (3/2)ACo\/££, that is, 
by apparently larger force for the compressive part of the 
impact (£ > 0 ) and apparently smaller for the restoring 
part (£ < 0 ). 


7 Conclusion 

We derive a new expression for the dissipative force acting 
between viscoelastic bodies during an impact. Contrary 
to the previous theories, based on the physically plausi¬ 
ble but non-rigorous approach, our theory exploits mathe¬ 
matically rigorous perturbation scheme with the small pa¬ 
rameter being the ratio of the microscopic relaxation time 
and the impact duration. We make calculations for the 
zero and first-order terms in this perturbation expansion. 
The new expression for the dissipative force noticeably 
differs from the previous one, obtained within the quasi¬ 
static approximation. Due to rigorous derivation from the 
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first principles we get a physically correct result for the 
dissipative force acting between bodies of different mate¬ 
rials; this was not possible within the previous approach. 
Moreover, our new theory is also lacking inconsistency of 
the previous theory with respect to materials with van¬ 
ishingly small elastic shear module. While the previous, 
quasi-static theory predicts the nonphysical zero dissipa¬ 
tion, the new theory implies dissipation, similar to that 
for ’’common” materials. 

In the present study we neglect the inertial effects, 
that is, we assume that the characteristic velocity of the 
problem is much smaller than the speed of sound in the 
bodies. The general approach presented in our study may 
be, however, further developed to take into account the 
inertial effects as well as high-order terms in the pertur¬ 
bation series. 

DSC and AVP acknowledge financial support from the Russian 
Scientific Foundation (grant no. 14-21-00090). 
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